// INCLUDES --------------------------------------------------------------------
#include "Sphere.h"

// FUNCTIONS -------------------------------------------------------------------

/* GENERAL:
   d    -> destination
   s    -> source
   size -> array size (0,1 for single elements)

   s    -> arbitrary sphere
   r    -> arbitrary ray
*/

// clear to c:[0,0,0,1], radius: 1.0f
void SphereClear(Sphere *s, ui32 size)
{
  GUARD_NULL("SphereClear *s", s);
  
  // clamp range to at least one element
  CLAMP_LOW(size, 1);

  // clear to c:[0,0,0,1], radius: 1.0f
  i32 i;
  for (i = 0; i < (i32)size; i++)
  {
    Vector4Clear(&s[i].c, 0);      // set c:[0,0,0,1]
    s[i].radius = 1.0f;            // set radius=1
    s[i].zNear = s[i].zFar = 0.0f; // reset ray hit distances to 0
  }
};            

// clone sphere src->dest
void SphereCopy(Sphere *d, Sphere *s, ui32 size)
{
  // copy contents
  FillCopy(d, s, sizeof(Sphere), size);
}; 


/* ray->sphere intersection testing

   a sphere with center c and radius r is described by:

   (x-cx)^2 + (y-cy)^2 + (z-cz)^2 - r^2 = 0

   every point that satisfies this equation is on the sphere.
   in vector form this becomes:
    
   (p - c)^2 - r^2 = 0
   
   if we now plug the ray-equation:
   
   p(t) = o + td

   as p into this, we get:

	 (o + td - c)^2 - r^2 = 0

	 which shows if the ray hits the sphere or not (means -> p(t) is on the sphere).
   rearranging this for t gives us:

	 (td + o - c)^2 - r^2 = 0

   now we expand the plugged-in ray equation via the binominal law:

   (a+b)^2 = a^2 + 2ab + b^2

   -> ((td)+(o-c))^2
      ((td)+(o-c))*((td)+(o-c))
      
      (td)(td)   = t^2 * d^2 
      (td)(o-c)  = tdo - tdc
      (o-c)(td)  = tdo - tdc
      (o-c)(o-c) = (o-c)^2   <- we leave this be for convenience

   ... compact the middle term a bit:

   2tdo - 2tdc -> 2d(o-c)t

   => d^2 * t^2 + 2d(o-c)t + (o-c)^2

   and combined with the rest of the formula this expands to the quadratic equasion:

	  a    x^2 +      p   x +        q          = 0

	 d^2 * t^2 + 2d(o - c)t + ((o - c)^2 - r^2)	= 0

   TODO: re-discover the re-arranging of the pq-formula for At^2+Bt+C = 0

	 which can be solved via the pq-formula for a standard quadratic equation,
	 giving us:

	        p           p^2 - 4aq
	 x1 = - -- + sqrt ( --------- )
	        2a            4a^2

	           and

          p           p^2 - 4aq
	 x2 = - -- - sqrt ( --------- )
	        2a            4a^2

	 eliminating the common factor of 2a in both fractions leaves us with:

	      -p + sqrt(p^2 - 4aq)
	 x1 = --------------------
	              2a

	           and

	      -p + sqrt(p^2 - 4aq)
	 x2 = --------------------
	              2a

	 all we have to do now is substitude all the variables with the right
	 algebraic terms:


	      -(2d(o - c)) + sqrt( (2d(o - c))^2 - 4 d^2 ((o - c)^2 - r^2) )
	 t =  ----------------------------------------------------------------
			  	                           2d


	                                  and


	      -(2d(o - c)) - sqrt( (2d(o - c))^2 - 4 d^2 ((o - c)^2 - r^2) )
	 t =  ----------------------------------------------------------------
			  	                           2d


	 the discriminant tells us if the sphere has been missed, grazed or
	 traversed.

	 D < 0 -> sphere has been missed (no solution in R)
	 D = 0 -> sphere has been grazed
   D > 0 -> sphere has been traversed by the ray, entering and leaving
            at the given u values
*/
// compute intersection
void SphereIntersect(Sphere *s, Ray *r, ui32 size)
{
  GUARD_NULL("SphereIntersect *s", s);
  GUARD_NULL("SphereIntersect *r", r);
  
  // clamp range to at least one element
  CLAMP_LOW(size, 1); 

  ui32 i;
  for (i = 0; i < size; i++)
  {
    Vector4 o_minus_c;
    Vector4Sub(&o_minus_c, &r->o, &s[i].c);                  // (o - c) = ray.origin - sphere.center

    fp a = Vector4ScalarProduct(&r->d, &r->d),               // a = d^2
       p = Vector4ScalarProduct(&r->d, &o_minus_c) * 2.0f,   // p = 2d(o - c)
       r = s[i].radius * s[i].radius,                        // r^2
       q = Vector4ScalarProduct(&o_minus_c, &o_minus_c) - r; // q = (o - c)^2 - r^2

	  fp discriminant = p*p - 4*a*q;                           // p^2 - 4aq

	  // check the discriminant for being >= 0
	  if (discriminant < 0) return;

    // ... and save the actual hit distances
	  discriminant = sqrt(discriminant);            // get the square root
   
    s[i].zNear = (-p - discriminant) / 2.0f*a;   // enter/graze distance

    if (discriminant != 0)
      s[i].zFar  = (-p + discriminant) / 2.0f*a; // exit distance
    else
      s[i].zFar = 0;                             // sphere grazed, no exit distance
  };
};
